Project 3 - Part 1-4

FFT and Aligning Images

import numpy as np
import matplotlib.pyplot as plt

# Returns the num channels for a passed in image np ndarray
def get_num_channels(img):
  # Check for grayscale
  if len(img.shape) < 3:
    return 1

  return img.shape[-1]

# Convert images to grayscale and print shape
def convert_grayscale(image):
  weights = [0.3, 0.6, 0.1]

  channels = get_num_channels(image)
  if channels == 1: # Grayscale
    return image
  elif channels == 3: # RGB
    return (image @ weights) / 255.0
  elif channels == 4: # ARGB
    # Only use the first 3 channels
    return (image[:,:,:3] @ weights) /255.0

Read in cell images and create other images from other image.

PATH = './cell_images'
files = [
    '0001.000.png',
    '0001.001.png',
    '0001.002.png',
    '0001.004.png',
    '0001.005.png',
    '0001.006.png',
]

images = []
for fname in files:
  image = plt.imread(f'{PATH}/{fname}')
  print(f'{fname}')
  print(f'shape: {image.shape}')
  print(f'max: {image.max()}')
  print(f'min: {image.min()}')
  print()

  # Convert images to grayscale
  images.append(convert_grayscale(image))

0001.000.png
shape: (510, 510)
max: 1.0
min: 0.0

0001.001.png
shape: (510, 510)
max: 1.0
min: 0.0

0001.002.png
shape: (510, 510)
max: 1.0
min: 0.0

0001.004.png
shape: (510, 510)
max: 1.0
min: 0.0

0001.005.png
shape: (510, 510)
max: 1.0
min: 0.0

0001.006.png
shape: (510, 510)
max: 1.0
min: 0.0

def display_images(imgs, fig_name):
  rows = int(np.ceil(len(images) / 4))
  fig, axs = plt.subplots(rows,4, figsize=(rows*4,4))
  fig.tight_layout()
  for i, axe in enumerate(axs):
    for j, ax in enumerate(axe):
      idx = (i*4) + j
      if idx < len(imgs):
        ax.imshow(imgs[idx], cmap='gray')
        ax.set_title(files[idx])
      # ax.set_xlabel(files[idx])

  plt.suptitle('Cell Images')
  plt.tight_layout()
  plt.savefig(fig_name)

display_images(images, 'gs_images.png')

Read in and create overlapping crops and images from the drone image as another example.

def create_crops_from_image(img_path):
    oth_images = []
    oimage = plt.imread(img_path)
    for i in range(10):
        x = np.random.randint(2000, 2400)
        y = np.random.randint(1000, 1400)
        
        oth_images.append(oimage[x:x+500, y:y+500])
        plt.imsave(f'./other_images/drone_00{i}.png', oth_images[i])
    
    # Convert images to grayscale to use
    gs_images = []
    for i in range(len(oth_images)):
        gs_images.append(convert_grayscale(oth_images[i]))
        
        if i < 3:
            plt.figure()
            plt.imshow(gs_images[i], cmap='gray')
        
    return gs_images
        
        
other_images = create_crops_from_image('./other_images/drone.jpg')

1) FFTs

I followed the book for how to do the low pass filter in the Fourier domain and display multiple images power spectrums, the filters, and the resulting image back in the spatial domain after it has been filtered in the Fourier domain. I implemented the main 3 types of low pass filters described in the textbook which are the ideal low pass filter (ILPF), gaussian low pass filter (GLPF), and Butterworth low pass filter (BLPF). The ideal low pass filter has a very sharp transition like a circle while the gaussian and Butterworth filters have additional parameters that can make a softer transition as shown in the examples below.

def ddistance(u,v,p,q):
  return (((u - (p/2)) ** 2) + ((v - (q/2)) ** 2)) ** (1/2)

# Ideal low pass filter of size (p,q)
def ideal_low_pass_filter(d0, p, q):
  out = np.zeros((p,q))

  for i in range(p):
    for j in range(q):
      out[i,j] = 1 if ddistance(i,j,p,q) <= d0 else 0

  return out

# Gaussian low pass filter of size (p,q)
def gaussian_low_pass_filter(d0, p, q):
  out = np.zeros((p,q))

  for i in range(p):
    for j in range(q):
      out[i,j] = np.exp((-1*(ddistance(i,j,p,q) ** 2)) / (2*((d0) ** 2)))

  return out

# Butterworth low pass filter of size (p,q)
def butterworth_low_pass_filter(d0, n, p, q):
  out = np.zeros((p,q))

  for i in range(p):
    for j in range(q):
      out[i,j] = (1) / (1 + (ddistance(i,j,p,q) / d0) ** (2*n))

  return out

Examples of the 3 low pass filters: Ideal, Gaussian, and Butterworth with different parameters so you can see what the filters look like. You can see the the Ideal LPF has a very sharp change while the Gaussian and Butterworth LPF can be more gradual depending on what parameters are passed to the filters.

def generate_low_pass_filters():
  
  ilpf_arr = [50, 100, 150]
  glpf_arr = [10, 50, 100]
  blpf_arr = [(50, 2), (50,4), (100,2)]
  
  rows = 3
  cols = 3
  fig, axs = plt.subplots(rows, cols, figsize=(10,8))
  for i, ax in enumerate(axs):
    
    ax[0].imshow(ideal_low_pass_filter(ilpf_arr[i], 1000, 1000), cmap='gray')
    ax[0].set_title(f'Ideal LPF: d0={ilpf_arr[i]}')
    
    ax[1].imshow(gaussian_low_pass_filter(glpf_arr[i], 1000, 1000), cmap='gray')
    ax[1].set_title(f'Gaussian LPF: d0={glpf_arr[i]}')
    
    ax[2].imshow(butterworth_low_pass_filter(blpf_arr[i][0], blpf_arr[i][1], 1000, 1000), cmap='gray')
    ax[2].set_title(f'Butterworth LPF: d0={blpf_arr[i][0]}, n={blpf_arr[i][1]}')


  plt.suptitle(f'Low Pass Filters')
  plt.tight_layout()
  plt.savefig(f'low-pass-filters.png')
  
generate_low_pass_filters()

Down below I display the original images, the power spectrum and the log of the power spectrum of the FFT of the images. I also display the filter and apply that filter in the Fourier domain. I do the padding because that is what is described in the course textbook in section (4.7 The Basics of Filtering in the Frequency Domain, see Figure 4.35 and above it where they describe the process of filtering in the Fourier Domain if needed). It sounds like it is not mandatory but it can be used to apply different effects to the image. You can see that I organize the frequencies so that zero is in the middle of the FFT. You don't see much only looking at the power spectrum, but you do see more information when you look at the log of the power spectrum which has some strong lines. I also show the Gaussian low pass filter that is used to filter in Fourier domain and the output image after the filtering. You can see that since it removes some of the higher frequencies that the image becomes slightly blurred and more soft.

# Multiply image by (-1)^(x+y) to center the Fourier transform
def center_fourier_transform(img_in):
  img = img_in.copy()
  x, y = img.shape
  for i in range(x):
    for j in range(y):
      img[i,j] = img[i,j] * ((-1) ** (i+j))

  return img


def pad_image(img):
  m,n = img.shape
  padded_img = np.zeros((2*m,2*n))
  padded_img[:m,:n] = img
  return padded_img


def display_power_spectrum(img, axs):
  # Display original image
  ax = axs[0]
  ax.imshow(img, cmap='gray')
  ax.set_title('Original')
  ax.set_xticks([])
  ax.set_yticks([])

  # Pad the image to size P and Q
  m,n = img.shape
  padded_img = pad_image(img)

  # Center and computer fourier transform
  fft = np.fft.fft2(center_fourier_transform(padded_img))

  # Magnitude of the DFT is the frequency spectrum
  # |F(u)| = [R(u)^2 + I(u)^2]^(1/2)
  ps = ((fft.real ** 2) + (fft.imag ** 2)) ** (1/2)

  # print(ps.shape)
  # print(ps.max())
  # print(ps.min())
  # print()

  # padded image
  ax = axs[1]
  ax.imshow(padded_img, cmap='gray')
  ax.set_title('Padded Image')
  ax.set_xticks([])
  ax.set_yticks([])

  # power spectrum, no log
  ax = axs[2]
  ax.imshow(ps, cmap='gray')
  ax.set_title('Power Spectrum')
  ax.set_xticks([])
  ax.set_yticks([])

  # power spectrum with log
  # s = c*log(1+r)
  ax = axs[3]
  ax.imshow(np.log(1 + ps), cmap='gray')
  ax.set_title('Log Power Spectrum')
  ax.set_xticks([])
  ax.set_yticks([])

  # gaussian low pass filter
  filter = gaussian_low_pass_filter(12, ps.shape[0], ps.shape[1])
  ax = axs[4]
  ax.imshow(filter, cmap='gray')
  ax.set_title('Gaussian Low Pass Filter')
  ax.set_xticks([])
  ax.set_yticks([])

  # Inverse FFT with filter
  out = center_fourier_transform(np.fft.ifft2(fft * filter).real)
  ax = axs[5]
  ax.imshow(out, cmap='gray')
  ax.set_title('Inverse FFT with Filter')
  ax.set_xticks([])
  ax.set_yticks([])

  # Out Image
  ax = axs[6]
  ax.imshow(out[:m, :n], cmap='gray')
  ax.set_title('Inverse FFT Image')
  ax.set_xticks([])
  ax.set_yticks([])

  return True

def generate_power_spectrum(imgs):
  rows = len(imgs)
  cols = 7
  fig, axs = plt.subplots(rows, cols, figsize=(16,14))
  for i, ax in enumerate(axs):
    idx = i
    img = imgs[idx]

    display_power_spectrum(img, ax)


  plt.suptitle(f'FFT Power Spectrum')
  plt.tight_layout()
  plt.savefig(f'fft-power-spectrum.png')

generate_power_spectrum(images)

Down below you can see the process as described above except on the other images. You can see that some of the log power spectrums for some of these images have additional vertical and horizontal lines for some of the stronger lines in the image. There is also a few other smaller lines in the log power spectrum that somewhat mirror other lines in the spatial domain. The Gaussian low pass filter has the same impact of removing some of the higher frequencies and blurring the image once we convert it back to the spatial domain.

generate_power_spectrum(other_images[:5])

The below chart shows the impact of the three different low pass filters in the Fourier domain with different parameters on different images. Its interesting to note the slight ringing effect the top ideal LPF has on the image because of the sharp edge in the filter. The lower Gaussian filter really blurs the image while the higher parameters for the Gaussian filter just slightly blurs the image. It's the same for the Butterworth filter, depending on what parameters you pass to the function.

def apply_filter(img, filter):
    # Pad the image to size P and Q
  m,n = img.shape
  padded_img = pad_image(img)

  # Center and computer fourier transform
  fft = np.fft.fft2(center_fourier_transform(padded_img))

  # Magnitude of the DFT is the frequency spectrum
  # |F(u)| = [R(u)^2 + I(u)^2]^(1/2)
  ps = ((fft.real ** 2) + (fft.imag ** 2)) ** (1/2)

  # Inverse FFT with filter
  out = center_fourier_transform(np.fft.ifft2(fft * filter).real)

  # Out Image
  return out[:m, :n]


def generate_low_pass_filters_on_image(img):
  
  ilpf_arr = [50, 100, 150]
  glpf_arr = [10, 50, 100]
  blpf_arr = [(50, 2), (50,4), (100,2)]
  
  rows = 3
  cols = 7
  fig, axs = plt.subplots(rows, cols, figsize=(16,8))
  for i, ax in enumerate(axs):
    
    ax[0].imshow(img, cmap='gray')
    ax[0].set_title('Original')
    
    m,n = img.shape
    p,q = m*2,n*2
    
    filter = ideal_low_pass_filter(ilpf_arr[i], p, q)
    ax[1].imshow(filter, cmap='gray')
    ax[1].set_title(f'ILPF: d0={ilpf_arr[i]}')
    
    ax[2].imshow(apply_filter(img, filter), cmap='gray')
    ax[2].set_title(f'ILPF: d0={ilpf_arr[i]}')
    
    filter = gaussian_low_pass_filter(glpf_arr[i], p, q)
    ax[3].imshow(filter, cmap='gray')
    ax[3].set_title(f'GLPF: d0={glpf_arr[i]}')
    
    ax[4].imshow(apply_filter(img, filter), cmap='gray')
    ax[4].set_title(f'GLPF: d0={glpf_arr[i]}')
    
    filter = butterworth_low_pass_filter(blpf_arr[i][0], blpf_arr[i][1], p, q)
    ax[5].imshow(filter, cmap='gray')
    ax[5].set_title(f'BLPF: d0={blpf_arr[i][0]}, n={blpf_arr[i][1]}')
    
    ax[6].imshow(apply_filter(img, filter), cmap='gray')
    ax[6].set_title(f'BLPF: d0={blpf_arr[i][0]}, n={blpf_arr[i][1]}')
    
    for j in range(cols):
      ax[j].set_xticks([])
      ax[j].set_yticks([])


  plt.suptitle(f'Low Pass Filters On Image')
  plt.tight_layout()
  plt.savefig(f'low-pass-filters-on-image.png')
  
generate_low_pass_filters_on_image(images[0])

generate_low_pass_filters_on_image(other_images[0])

2) Phase Correlation

I implemented the phase correlation with the low pass filter. I show examples of the two original images, the phase correlation output, and what it looks like with different low pass filters and parameters. You can see in the raw phase correlation that there is usually a pixel that is brighter than the rest that shows the offset between the images. You can see that this offset actually becomes easier to see with some of the Ideal, Gaussian, and Butterworth low pass filters. This is very obvious for the last three where the bright pixel becomes more obvious to see with the LPF. On the top two, you can see the brightest pixel on the top left and then on the second one the brightest pixel is on the bottom right.

# Get the index of the max value
def get_max_index(pc):
  max_idx = None
  for i in range(pc.shape[0]):
    for j in range(pc.shape[1]):
      if max_idx is None:
        max_idx = (i,j)
        max_val = pc[i,j]
      elif pc[i, j] >= max_val:
        max_idx = (i,j)
        max_val = pc[i,j]

  return max_idx, max_val

def phase_correlation(f, g, ftype=0, filter_params={}):
  m,n = f.shape
  
  # Center for low pass filter
  f_fft = np.fft.fft2(center_fourier_transform(f))
  g_fft = np.fft.fft2(center_fourier_transform(g))

  # G* = R(u) - jI(u)
  g_star_fft = g_fft.real - g_fft.imag
  a = f_fft * g_star_fft

  # |F(u)|=[(R(u) ** 2) + (I(u) ** 2)] ** (1/2)
  mag = np.sqrt((a.real ** 2) + (a.imag ** 2))
  b = a/mag

  # Use low pass filter
  if ftype == 1:
    b = b * ideal_low_pass_filter(filter_params['d0'], m, n)
  elif ftype == 2:
    b = b * gaussian_low_pass_filter(filter_params['d0'], m, n)
  elif ftype == 3:
    b = b * butterworth_low_pass_filter(filter_params['d0'], filter_params['n'], m, n)

  return center_fourier_transform(np.fft.ifft2(b).real)

def generate_phase_correlation(img_pairs):
  rows = len(img_pairs)
  cols = 6
  fig, axs = plt.subplots(rows, cols, figsize=(16,14))
  for i, ax in enumerate(axs):
    idx = i
    f = img_pairs[idx][0]
    g = img_pairs[idx][1]

    ax[0].imshow(f, cmap='gray')
    ax[0].set_title('f image')
    ax[1].imshow(g, cmap='gray')
    ax[1].set_title('g image')
    
    ax[2].imshow(phase_correlation(f,g), cmap='gray')
    ax[2].set_title('PC - no filter')
    
    ax[3].imshow(phase_correlation(f,g,2,{'d0':50}), cmap='gray')
    ax[3].set_title(f'PC - ILPF - d0=50')
    
    ax[4].imshow(phase_correlation(f,g,2,{'d0':50}), cmap='gray')
    ax[4].set_title(f'PC - GLPF - d0=50')
    
    ax[5].imshow(phase_correlation(f,g,3,{'d0':50,'n':3}), cmap='gray')
    ax[5].set_title(f'PC - BLPF - d0=50,n=3')
    
    for j in range(cols):
      ax[j].set_xticks([])
      ax[j].set_yticks([])


  plt.suptitle(f'Phase Correlation')
  plt.tight_layout()
  plt.savefig(f'phase-correlation.png')
  
img_pairs = [(images[4], images[5]), (images[1], images[0]), (other_images[0], other_images[1]), (other_images[2], other_images[3]), (other_images[0], other_images[3])]
generate_phase_correlation(img_pairs)

Down below I display a comparison on changing the parameters to the filters when doing the phase correlation and the output. You can see that as you change the parameters to the filters is really does highlight where the that peak in the phase correlation is, so there is probably an optimal low pass params for certain types of images when doing the filtering with the phase correlation.

def param_to_string(param):
  # {'d0': 10, 'n': 2}
  s = f'd0={param["d0"]}'
  if param.get('n', '') != '':
    s += f',n={param["n"]}'
    
  return s

def generate_phase_correlation_params(f,g):
  
  fnames = ['ILPF', 'GLPF', 'BLPF', 'BLPF_1', 'BLPF_2']
  ftypes = [1, 2, 3, 3, 3]
  fparams = {
    'ILPF': [{'d0': 25},{'d0': 100},{'d0': 200}],
    'GLPF': [{'d0': 25},{'d0': 50},{'d0': 100}],
    'BLPF': [{'d0': 25, 'n': 1},{'d0': 50, 'n': 1},{'d0': 100, 'n': 1}],
    'BLPF_1': [{'d0': 25, 'n': 3},{'d0': 50, 'n': 3},{'d0': 100, 'n': 3}],
    'BLPF_2': [{'d0': 25, 'n': 5},{'d0': 50, 'n': 5},{'d0': 100, 'n': 5}],
  }
  
  rows = len(fnames)
  cols = 6
  fig, axs = plt.subplots(rows, cols, figsize=(16,14))
  for i, ax in enumerate(axs):
    idx = i

    ax[0].imshow(f, cmap='gray')
    ax[0].set_title('f image')
    ax[1].imshow(g, cmap='gray')
    ax[1].set_title('g image')
    
    fname = fnames[i]
    ftype = ftypes[i]
    
    ax[2].imshow(phase_correlation(f,g), cmap='gray')
    ax[2].set_title('PC - no filter')
    
    fparam = fparams[fname][0]
    ax[3].imshow(phase_correlation(f,g,ftype,fparam), cmap='gray')
    ax[3].set_title(f'PC - {fname} - {param_to_string(fparam)}')
    
    fparam = fparams[fname][1]
    ax[4].imshow(phase_correlation(f,g,ftype,fparam), cmap='gray')
    ax[4].set_title(f'PC - {fname} - {param_to_string(fparam)}')
    
    fparam = fparams[fname][2]
    ax[5].imshow(phase_correlation(f,g,ftype,fparam), cmap='gray')
    ax[5].set_title(f'PC - {fname} - {param_to_string(fparam)}')
    
    for j in range(cols):
      ax[j].set_xticks([])
      ax[j].set_yticks([])


  plt.suptitle(f'Phase Correlation - Filter Params')
  plt.tight_layout()
  plt.savefig(f'phase-correlation-filter-params.png')


generate_phase_correlation_params(other_images[0], other_images[1])

3) Peak Finding

I wrote out the functionality to perform the phase correlation and then find the pixel location with the maximal value. To then illustrate whether that implementation was correct I also align the images by resolving the four cases of ambiuity (as described in part 4 Mosaic Building). That way I can align the image using the peak found in the phase correlation and display the images together. You can see from the examples below that for my images using the pixel location with the maximal value works extremely well for aligning these images and it is very satisfying to see the images properly align. In my testing for my images, it seemed like using a simple threshold of 0.01 worked best for finding the peak value of the phase correlation and determining whether or not the images overlap. In practice, looking at the connected components from the phase correlation image thresholded at some value and then averaging the magnitude and position within those components in order to find the best peak to use for the two images works the best if you have correct thresholds. This is especially true depending on what low pass filter you use during the phase correlation and what parameters you pass to that filter. For example, if you look at the examples above you can see that certain filters make it so that the bright area in the phase correlation is a lot larger and consists of multiple pixels. In this scenario the base thing to do would be to use the thresholded connected components and do averaging on the position within those components find the best peak within that patch. Overall for most scenarios it seems like just using the pixel location with the maximal value performs very well as a simpler option to tune and use for these images.

import skimage

def pad_image_with_average(f,g):
  f_temp = f.copy()
  g_temp = g.copy()
  
  fm, fn = f_temp.shape
  gm, gn = g_temp.shape
  if fm > gm:
    temp = np.ones((fm,gn)) * g.mean()
    temp[:gm, :gn] = g_temp
    g_temp = temp
  
  fm, fn = f_temp.shape
  gm, gn = g_temp.shape
  if gm > fm:
    temp = np.ones((gm,fn)) * f.mean()
    temp[:fm, :fn] = f_temp
    f_temp = temp
  
  fm, fn = f_temp.shape
  gm, gn = g_temp.shape
  if fn > gn:
    temp = np.ones((gm,fn)) * g.mean()
    temp[:gm, :gn] = g_temp
    g_temp = temp
  
  fm, fn = f_temp.shape
  gm, gn = g_temp.shape
  if gn > fn:
    temp = np.ones((fm,gn)) * f.mean()
    temp[:fm, :fn] = f_temp
    f_temp = temp
    
  return f_temp, g_temp


def normalized_correlation(a,b):
  a = a - a.mean()
  b = b - b.mean()

  a_m = np.sqrt(np.square(a).sum())
  b_m = np.sqrt(np.square(b).sum())

  return ((a*b)/(a_m*b_m)).sum()


def find_best_corner(f, g, max_idx, display=True):
  m, n = f.shape
  gm, gn = g.shape
  best = None

  # Setup second image
  out2 = np.zeros((m*2,n*2))
  out2[max_idx[0]:max_idx[0]+gm, max_idx[1]:max_idx[1]+gn] = g

  # Try the best
  out1 = np.zeros((m*2,n*2))
  out1[:m,:n] = f
  out1[m:2*m,:n] = f
  out1[:m,n:2*n] = f
  out1[m:2*m,n:2*n] = f

  tl1 = out1[max_idx[0]:m, max_idx[1]:n]
  tr1 = out1[m:max_idx[0]+m, max_idx[1]:n]
  bl1 = out1[max_idx[0]:m, n:max_idx[1]+n]
  br1 = out1[m:max_idx[0]+m, n:max_idx[1]+n]
  
  tl2 = out2[max_idx[0]:m, max_idx[1]:n]
  tr2 = out2[m:max_idx[0]+m, max_idx[1]:n]
  bl2 = out2[max_idx[0]:m, n:max_idx[1]+n]
  br2 = out2[m:max_idx[0]+m, n:max_idx[1]+n]

  tl = normalized_correlation(tl1,tl2)
  tr = normalized_correlation(tr1,tr2)
  bl = normalized_correlation(bl1,bl2)
  br = normalized_correlation(br1,br2)

  if display:
    print('\n4 Patches Correlations:')
    print(tl)
    print(tr)
    print(bl)
    print(br)

  arr = [tl, tr, bl, br]
  best_idx = np.argmax(arr)
  
  out = np.zeros((m*2,n*2))
  if best_idx == 0:
    out[:m,:n] = f
  elif best_idx == 1:
    out[m:2*m,:n] = f
  elif best_idx == 2:
    out[:m,n:2*n] = f
  elif best_idx == 3:
    out[m:2*m,n:2*n] = f


  out[max_idx[0]:max_idx[0]+gm, max_idx[1]:max_idx[1]+gn] = g
  best = out

  return best, arr[best_idx], best_idx


def get_fitted_image(inp_img):
    mask = (inp_img > 0).astype(int)
    min_i = np.inf
    min_j = np.inf
    max_i = -1 * np.inf
    max_j = -1 * np.inf
    for i in range(inp_img.shape[0]):
        for j in range(inp_img.shape[1]):
            if mask[i,j] > 0:
                if i < min_i:
                    min_i = i
                if j < min_j:
                    min_j = j
            if mask[i,j] > 0:
                if i > max_i:
                    max_i = i
                if j > max_j:
                    max_j = j

    fitted_canvas = inp_img[min_i:max_i,min_j:max_j]
    return fitted_canvas


def image_registration(f, g, display=True):
  f_temp, g_temp = pad_image_with_average(f, g)
  
  if display:
    fig, axs = plt.subplots(1, 5, figsize=(14,4))
    axs[0].imshow(f, cmap='gray')
    axs[0].set_title('f image')
    
    axs[1].imshow(g, cmap='gray')
    axs[1].set_title('g image')
  
  # Compute the phase correlation
  pc = phase_correlation(f_temp, g_temp)
  
  # Simple peak finding - max value index
  max_idx, max_val = get_max_index(pc)
  if display:
    print(f'max_val:{max_val}')
    if max_val < 0.01:
      print(f'max val is less than threshold of 0.01, images dont overlap')
  
  # Display phase correlation and max peak
  if display:
    # print(f'max: ({max_idx[0]},{max_idx[1]})')
    axs[2].imshow(pc, cmap='gray')
    axs[2].set_title('phase correlation with argmax')
    # imshow flips an axis and plot doesnt
    axs[2].plot(max_idx[1], max_idx[0], 'ro')

  # Resolve the four cases of ambiguity after you find the peek
  # based on the correlation on the overlapping patches in the
  # spatial domain as described in Mosaic Building
  best, best_corr, corner = find_best_corner(f, g, max_idx, display)
  if display:
    axs[3].imshow(best, cmap='gray')
    axs[3].set_title('Aligned Images')
  
  fitted_best = get_fitted_image(best)
  if display:
    axs[4].imshow(fitted_best, cmap='gray')
    axs[4].set_title('Fitted Aligned Images')

  if display:
    for i in range(5):
      axs[i].set_xticks([])
      axs[i].set_yticks([])
      
    plt.suptitle(f'Peak Finding - Argmax')
    # plt.tight_layout()
    plt.savefig(f'peak-finding-argmax.png')

  return best, best_corr, max_idx, corner


def connected_components_peak_finding(pc, threshold=0.01):
  print(pc.max())
  pc_binary = (pc > threshold).astype(int)
  print(f'num pixels greater than {threshold}: {pc_binary.sum()}')
  
  # Connected components using the thresholded image
  labels, num_labels = skimage.measure.label(pc_binary, 0, True)
  print(f'num labels: {num_labels}')
  
  # Get the region props
  props = skimage.measure.regionprops(labels, pc)
  print(f'len of props: {len(props)}')
  
  # Round to get correct index
  rounded = np.round(props[0].centroid_weighted)
  
  # Using centroid weighted in order to get the centroid weighted with the intensity range.
  return (int(rounded[0]),int(rounded[1]))


def image_registration_connected_components(f, g, display=True):
  f_temp, g_temp = pad_image_with_average(f, g)
  
  if display:
    fig, axs = plt.subplots(1, 5, figsize=(14,4))
    axs[0].imshow(f, cmap='gray')
    axs[0].set_title('f image')
    
    axs[1].imshow(g, cmap='gray')
    axs[1].set_title('g image')
  
  # Compute the phase correlation
  pc = phase_correlation(f_temp, g_temp, 1, {'d0': 200})
      
  # Connected components peak finding
  # if not argmax_peak_finding:
  max_idx = connected_components_peak_finding(pc)
  max_val = pc[max_idx[0],max_idx[1]]
  print(max_idx)
  print(max_val)
  if max_val < 0.01:
    print(f'max val is less than threshold of 0.01, images dont overlap')
  
  # Display phase correlation and max peak
  if display:
    # print(f'max: ({max_idx[0]},{max_idx[1]})')
    axs[2].imshow(pc, cmap='gray')
    axs[2].set_title('phase correlation with argmax')
    # imshow flips an axis and plot doesnt
    axs[2].plot(max_idx[1], max_idx[0], 'ro')

  # Resolve the four cases of ambiguity after you find the peek
  # based on the correlation on the overlapping patches in the
  # spatial domain as described in Mosaic Building
  best, best_corr, corner = find_best_corner(f, g, max_idx, display)
  if display:
    axs[3].imshow(best, cmap='gray')
    axs[3].set_title('Aligned Images')
  
  fitted_best = get_fitted_image(best)
  if display:
    axs[4].imshow(fitted_best, cmap='gray')
    axs[4].set_title('Fitted Aligned Images')

  if display:
    for i in range(5):
      axs[i].set_xticks([])
      axs[i].set_yticks([])
      
    plt.suptitle(f'Peak Finding - Connected Components')
    # plt.tight_layout()
    plt.savefig(f'peak-finding-connected-components.png')

  return best, best_corr, max_idx, corner

Examples of peak finding by using the pixel location with the maximal value. I also printed and tested different thresholds for determining if the two images actually overlap. These examples below are specifically images that do overlap. I found that a value of 0.01 seemed to be a good threshold for highest pixel value for determining whether or not the images actually overlap. You can see the max values below which are between 0.02 - 0.26.

out = image_registration(images[4], images[5])

max_val:0.026673842529438815

4 Patches Correlations:
5.9705083025726226e-05
0.966170162052822
-0.1658247545475127
0.3542002832093809

out = image_registration(images[1], images[0])

max_val:0.05321337484149325

4 Patches Correlations:
0.0458111236335462
0.08873192701294352
0.962104436797337
0.07025478570448268

out = image_registration(other_images[1], other_images[4])

max_val:0.2252031832125647

4 Patches Correlations:
0.17668686030474332
0.03198615584644228
-0.05364424538843698
1.0

out = image_registration(images[1],images[4])

max_val:0.05893441811363937

4 Patches Correlations:
0.9684790749779867
0.3019129044452967
-0.05244891206039295
-0.06227593086926022

NOTE: Examples of Images that do not overlap below.

Examples of images that do not overlap for finding a good threshold of maximum value to be below a certain threshold. You can see that for these three examples that the max value is only up to around 0.007 which is under the threshold of 0.01. NOTE: I still display the images aligned improperly just to demonstrate that they do not overlap properly, but I print out that the images do not overlap based on that max val being below the threshold.

out = image_registration(images[1],other_images[3])

max_val:0.006266453972425221
max val is less than threshold of 0.01, images dont overlap

4 Patches Correlations:
-0.006340864909351679
-0.029963400936686216
0.003736807406686117
-0.05463664066962649

out = image_registration(images[0],other_images[0])

max_val:0.0065026941749681145
max val is less than threshold of 0.01, images dont overlap

4 Patches Correlations:
0.03840708352365752
0.11247993084037902
-0.06313805591671455
0.0478074844440424

out = image_registration(images[4],other_images[2])

max_val:0.006687968061909436
max val is less than threshold of 0.01, images dont overlap

4 Patches Correlations:
0.11425482350366958
-0.0066857402146745515
-0.08631308961056362
0.054744001984065496

Examples of peak finding using connected components above a threshold of 0.01 and then averaging the magnitude and position within the found components in order to find the best peak for the offset between the two images.

out = image_registration_connected_components(images[4], images[5])

0.017702452726660628
num pixels greater than 0.01: 3
num labels: 1
len of props: 1
(107, 3)
0.017702452726660628

4 Patches Correlations:
5.9705083025726226e-05
0.966170162052822
-0.1658247545475127
0.3542002832093809

out = image_registration_connected_components(images[1],images[4])

0.04114481544165318
num pixels greater than 0.01: 5
num labels: 1
len of props: 1
(5, 360)
0.04114481544165318

4 Patches Correlations:
0.9684790749779867
0.3019129044452967
-0.05244891206039295
-0.06227593086926022

4) Mosaic Building

I created a function called mosaic_image_from_file() which will read in a list of images from a file, and then will create a mosaic from them and output them to a file. I also display the mosaic'd image in the notebook. The mosaic_image_from_file uses the mosaic_images function which computes the matrix of phase correlation and cross correlation from the aligned images and then uses that in order to align and put all images onto the canvas. You can see the mosaic created from the cell images and the beach trees images down below. In order to align the images properly, I use the technique in the last paragraph of (4 - Mosaic Building) on the assignment and resolve the four cases of ambiguity after the peak in phase correlation is found. I then computer the correlation on the overlapping patches in the spatial domain as described in the image on the homework. This helps resolve the four cases of ambiguity and figure out how to align the images.

def mosaic_images(imgs):
    # compute matrix of cross correlation values
    curr_images = imgs
    out_imgs = []
    corrs = np.zeros((len(curr_images), len(curr_images)))
    offset_arr = []
    corners = np.zeros((len(curr_images), len(curr_images)))
    for i in range(len(curr_images)):
        for j in range(len(curr_images)):
            out, corr, offset, corner = image_registration(curr_images[i], curr_images[j], False)
            out_imgs.append(out)
            corrs[i,j] = corr
            offset_arr.append(offset)
            corners[i,j] = corner
    
    print('matrix of cross correlation values for image pairs')
    print(corrs)
    
    # choose anchor image
    anchor_idx = np.argmax(corrs.sum(axis=1))
    print(f'anchor image idx: {anchor_idx}')
    
    # put images on canvas using anchor image as base
    m,n = curr_images[anchor_idx].shape
    canvas = np.zeros((m*3,n*3))
    so = m
    canvas[so:so+m,so:so+n] = curr_images[anchor_idx]
    for i in range(len(curr_images)):
        if i != anchor_idx:
            corner = corners[anchor_idx,i]
            offset = offset_arr[(corners.shape[1]*anchor_idx)+i]
            
            if corner == 0:
                canvas[so+offset[0]:so+m+offset[0],so+offset[1]:so+n+offset[1]] = curr_images[i]
            elif corner == 1:
                canvas[so+offset[0]-m:so+m+offset[0]-m,so+offset[1]:so+n+offset[1]] = curr_images[i]
            elif corner == 2:
                canvas[so+offset[0]:so+m+offset[0],so+offset[1]-n:so+n+offset[1]-n] = curr_images[i]
            elif corner == 3:
                canvas[so+offset[0]-m:so+m+offset[0]-m,so+offset[1]-n:so+n+offset[1]-n] = curr_images[i]

    fitted_canvas = get_fitted_image(canvas)
    
    fig, axs = plt.subplots(1,2, figsize=(8,6))
    axs[0].imshow(canvas, cmap='gray')
    axs[0].set_title('Mosaic')
    axs[1].imshow(fitted_canvas, cmap='gray')
    axs[1].set_title('Fitted Mosaic')
    
    plt.suptitle(f'Mosaic Building')
    # plt.tight_layout()
    plt.savefig(f'mosaic-building.png')
    
    return canvas, fitted_canvas


def mosaic_image_from_file(path, fname, out_file):
    imgs = []
    with open(fname, 'r') as f:
        for l in f:
            ipath = path + l.strip()
            imgs.append(convert_grayscale(plt.imread(ipath)))
            
    mosaic, fitted_mosaic = mosaic_images(imgs)
    plt.imsave(path + out_file, fitted_mosaic)
    return
            

Below is an example of building the mosaic from images already read in an array in memory. You can see that it computes the cross correlation matrix and then uses that to select the anchor image and build the mosaic from there. I resolve the four cases of ambiguity as described above and in the assignment information.

out = mosaic_images(images)

matrix of cross correlation values for image pairs
[[1.         0.96210444 0.06512203 0.97111869 0.97253864 0.05555162]
 [0.96210444 1.         0.95245673 0.97628675 0.96847907 0.96525974]
 [0.02836406 0.95245673 1.         0.10822695 0.98047127 0.92242006]
 [0.97111869 0.97628675 0.21037403 1.         0.94213967 0.0608977 ]
 [0.97253864 0.96847907 0.98047127 0.94213967 1.         0.96617016]
 [0.1322312  0.96525974 0.92242006 0.08987814 0.96617016 1.        ]]
anchor image idx: 4

out = mosaic_images(other_images)

C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:39: RuntimeWarning: Mean of empty slice.
  a = a - a.mean()
C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:40: RuntimeWarning: Mean of empty slice.
  b = b - b.mean()

matrix of cross correlation values for image pairs
[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
anchor image idx: 0

Below is the example of building a function that reads in a set of images specified in a file and builds and then saves a mosaic out to the desired location, based on the parameters that were passed into the function.

mosaic_image_from_file('./cell_images/', './cell_images/read_cells.txt', 'cell_mosaic.png')

C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:39: RuntimeWarning: Mean of empty slice.
  a = a - a.mean()
C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:40: RuntimeWarning: Mean of empty slice.
  b = b - b.mean()

matrix of cross correlation values for image pairs
[[1.         0.96210444 0.06512203 0.97111869 0.97253864 0.05555162]
 [0.96210444 1.         0.95245673 0.97628675 0.96847907 0.96525974]
 [0.02836406 0.95245673 1.         0.10822695 0.98047127 0.92242006]
 [0.97111869 0.97628675 0.21037403 1.         0.94213967 0.0608977 ]
 [0.97253864 0.96847907 0.98047127 0.94213967 1.         0.96617016]
 [0.1322312  0.96525974 0.92242006 0.08987814 0.96617016 1.        ]]
anchor image idx: 4

Below is the example of building a function that reads in a set of images specified in a file for the drone image and builds and then saves a mosaic out to the desired location, based on the parameters that were passed into the function.

mosaic_image_from_file('./other_images/', './other_images/drone.txt', 'drone_mosaic.png')

C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:39: RuntimeWarning: Mean of empty slice.
  a = a - a.mean()
C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:40: RuntimeWarning: Mean of empty slice.
  b = b - b.mean()

matrix of cross correlation values for image pairs
[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
anchor image idx: 0